splev Subroutine

      pure subroutine splev(t,n,c,k,x,y,m,e,ier)

      !  calling sequence:
      !     call splev(t,n,c,k,x,y,m,e,ier)
      !
      !  input parameters:
      !    t    : array,length n, which contains the position of the knots.
      !    n    : integer, giving the total number of knots of s(x).
      !    c    : array,length n, which contains the b-spline coefficients.
      !    k    : integer, giving the degree of s(x).
      !    x    : array,length m, which contains the points where s(x) must be evaluated.
      !    m    : integer, giving the number of points where s(x) must be evaluated.
      !    e    : integer, boundary condition for points outside the support
      !           0 = the spline is extrapolated from the end spans
      !           1 = the spline evaluates to zero for those points,
      !           2 = extrapolation not allowed, ier is set to 1 and the subroutine returns,
      !           3 = the spline evaluates to the value of the nearest boundary point.
      !
      !  output parameter:
      !    y    : array,length m, giving the value of s(x) at the different points.
      !    ier  : error flag
      !
      !  restrictions:
      !    m >= 1
      !
      !  other subroutines required: fpbspl.
      !
      !  references :
      !    de boor c  : on calculating with b-splines, j. approximation theory 6 (1972) 50-62.
      !    cox m.g.   : the numerical evaluation of b-splines, j. inst. maths applics 10 (1972) 134-149.
      !    dierckx p. : curve and surface fitting with splines, monographs on numerical analysis, oxford
      !                 university press, 1993.
      !
      !  author :
      !    p.dierckx
      !    dept. computer science, k.u.leuven
      !    celestijnenlaan 200a, b-3001 heverlee, belgium.
      !    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
      !
      !  ..scalar arguments..
      integer,     intent(in)  :: n, k, m, e
      integer,     intent(out) :: ier
      !  ..array arguments..
      real(RKIND), intent(in)  :: t(n), c(n), x(m)
      real(RKIND), intent(out) :: y(m)

      !  ..local scalars..
      integer :: i, k1, l, l1, nk1,k2
      real(RKIND) :: arg, tb, te
      !  ..local array..
      real(RKIND) :: h(MAX_ORDER+1)
      !  ..
      !  before starting computations a data check is made. if the input data
      !  are invalid control is immediately repassed to the calling program.
      ier = FITPACK_INPUT_ERROR
      if (m<1) return

      ier = FITPACK_OK
      !  fetch tb and te, the boundaries of the approximation interval.
      k1  = k  + 1
      k2  = k1 + 1
      nk1 = n - k1
      tb  = t(k1)
      te  = t(nk1 + 1)
      l   = k1
      l1  = l + 1

      !  main loop for the different points.
      user_points: do i = 1, m

        ! fetch a new x-value arg.
        arg = x(i)

        ! check if arg is in the support
        outside: if (arg<tb .or. arg>te) then
            select case (e)
               case (OUTSIDE_EXTRAPOLATE)
                ! Continue normally
               case (OUTSIDE_ZERO)
                  y(i) = zero
                  cycle user_points
               case (OUTSIDE_NOT_ALLOWED)
                  ier = FITPACK_INVALID_RANGE
                  return
               case (OUTSIDE_NEAREST_BND)
                  arg = max(min(arg,te),tb)
            end select
        endif outside

        ! search for knot interval t(l) <= arg < t(l+1)
        do while (arg<t(l) .and. l1/=k2)
          l1 = l
          l  = l - 1
        end do
        do while (arg>=t(l1) .and. l/=nk1)
          l  = l1
          l1 = l + 1
        end do

        ! evaluate the non-zero b-splines at arg.
        h = fpbspl(t, n, k, arg, l)

        ! find the value of s(x) at x=arg.
        y(i) = dot_product(c(l-k:l),h(1:k1))
      end do user_points

      end subroutine splev